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Abstract 

We are developing a framework for multiscale computation which 
enables models at a "microscopic" level of description, for example 
Lattice Boltzmann, Monte Carlo or Molecular Dynamics simulators, 
to perform modelling tasks at "macroscopic" length scales of interest. 
The plan is to use the microscopic rules restricted to small "patches" 
of the domain, the "teeth", using interpolation to bridge the "gaps". 
Here we explore general boundary conditions coupling the widely sep- 
arated "teeth" of the microscopic simulation that achieve high order 
accuracy over the macroscale. We present the simplest case when the 
microscopic simulator is the quintessential example of a partial differ- 
ential equation. We argue that classic high-order interpolation of the 
macroscopic field provides the correct forcing in whatever boundary 
condition is required by the microsimulator. Such interpolation leads 
to Tooth Boundary Conditions which achieve arbitrarily high-order 
consistency. The high-order consistency is demonstrated on a class of 
linear partial differential equations in two ways: firstly through the 
eigenvalues of the scheme for selected numerical problems; and sec- 
ondly using the dynamical systems approach of holistic discretisation 
on a general class of linear pdes. Analytic modelling shows that, for a 
wide class of microscopic systems, the subgrid fields and the effective 
macroscopic model are largely independent of the tooth size and the 
particular tooth boundary conditions. When applied to patches of 
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microscopic simulations these tooth boundary conditions promise ef- 
ficient macroscale simulation. We expect the same approach will also 
accurately couple patch simulations in higher spatial dimensions. 

Keywords: multiscale computation, gap tooth scheme, coupling boundary 
conditions, high order consistency 
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1 Introduction 

The components of physical systems often o perate on yastly different space 
and time scales fiDolbow. Khaleel fc Mitchell 2 004) . We must somehow sim- 
ulate such systems on the scale of interest and operation. But systems that 
depend on physical processes at multiple scales pose notorious difficulties. 
These multiscale difficulties are major obstacles to progress in fields as di- 
yerse as enyironmental and geosciences, climate, materials, combustion, high 
energy density ph ysics, fusion, bioscie nce, chemistry, power grids and infor- 



mation networks (jPolbow et al.!!2004! ) 



He re we further develop the equation free approach to multiscale mod- 
elling (jKeyrekidis et al l!2003!) . Giyen a numerical simulator for physical com- 



ponents at much smaller scales than the scale of primary interest, the aim 
of the methodology is to bridge the space and time scales to simulations re- 
solving the macroscale of interest. Here we focus on bridging space scales 
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by i mproving the accuracy of the gap-too t h met hodology for microsimula- 
tors (|Gear et all 1200.1 ISamaev et allboni 1200,^ . Crucially, our gap-tooth 
methods must adapt to whatever microsimulator code is provided; one key 
application of this work is to microsimulators that are tried and tested legacy 
codes that we do not want to modify. 

The equation-free approach provides on the fly closure methods which 
const itute critical components of, for example, mathematical homogeniza- 
tion ( Samaey et al. 200,4 Gustafsson fc Mossinol 2003LjBalakotaiah fc Chang 
'2003, e.g.), renormalization group techniques (jEi et alboOOLlMudavanhu 0'Ma11evl 
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Chorin Stinis"200.4 e.g.), and multiscale finite elements (iHc 
Chen fc Hou.2002. e.g.). These closure methods not only need to be 



computationally efficient but also need to be capable of reproducing the phys- 
ical dynamics with high fidelity. That is, we seek a methodology that can be 
systematically refined. 

Using microscopic simulators of the one dimensional Burgers' equation, 
Roberts fc KevrekidisI ^00^ demonstrated the possibility of achieving high 
order accuracy in the gap-tooth scheme for macroscale dynamics. The par- 
ticular microsimulator we use is a fine scale discretization of the pde which 
we execute only in the interior of the teeth (see Figure Q). At each time step 
during execution, the microsimulator within each tooth requires boundary 
values which must be continuously updated. If the microsimulator was to be 
executed over the entire macrodomain, these boundary values would natu- 
rally come from the immediately neighboring fine grid; this grid is missing 
in gap-tooth simulation. That pilot study only considered microsimulators 
which had boundary conditions of specified flux at the edges of their simula- 
tion teeth. Here we generalise the analysis to consider microsimulators with 
either 

• Dirichlet boundary conditions of specified field u at the tooth edges. 
Section |21 

• mixed boundary conditions of specified avj ±b9xVj at the tooth edges. 
Section ini or 

• nonlocal two-point boundary conditions such as those arising in a mi- 
croscale discretisation of a pde. Section |31 



Consider the gap-tooth scheme ( Gear et al. 2003L Samaev et al. 2004 . 
e.g.) illustrated in Figure [U Let Vj(x, t) be the fine scale, microscopic field 
in the jth tooth, and Uj the jth coarse grid value; that is, the value at the 
center of each tooth. Let the tooth width be h. Then the edge of a tooth 
lies at a distance h/2 from its coarse grid point, a fraction r = h./(2H] to 
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Figure 1: simulation of Burgers' equation using Dirichlet boundary condi- 
tions on the teeth (specified u on the edges). 

the next coarse grid point. The amount of computation performed by mi- 
crosimulators is proportional to the width of the (microscale) teeth. Hence 
we aim for the fraction r to be as small as possible, so that the teeth are 
a relatively small part of the physical domain and the computational cost 
minimised. The coupling rule developed in Sections El El and ID is that you 
obtain whatever values are necessary for the boundaries of the microscopic 
simulators by classic interpolation of the macroscopic grid values from neigh- 
bouring teeth. As a nonlinear example of the coupling we develop, Figure ^ 
shows a gap-tooth simulation of the nonlinear dynamics of Burgers' equation 
in one spatial dimension. 

This coupling rule promotes a strong connection between classic finite dif- 
ference discretisations of pdes, classic finite elements, and the methodology 
of the gap-tooth scheme. First, the pde acts as the quintessential example 
of a microsimulator in that it informs us of the dynamics in an 'infinitesimal 
patch'. The only difference between the gap-tooth scheme and the spatial 
discretisation of pdes is that the microsimulators in the gap-tooth scheme 
encode the dynamics on small finite patches, whereas the pde encodes the 
dynamics on infinitesimal patches. Consequently, classic interpolation serves 
the same role in both: namely, the interpolation appropriately transfers infor- 
mation from the macroscale of interest to the microscale simulators. Second, 
the theoretical support for the gap-tooth scheme is based upon a subgrid scale 
structure, as is the classic finite element method. Also, the solvability condi- 
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tion in the construction of the theoretical gap-tooth model is similar to the 
Galerkin projection of finite elements. But whereas finite elements impose a 
class of subgrid fields, both the theoretical approach here and the gap-tooth 
scheme use actual subgrid scale dynamics, obtained from the microsimulator 
or the PDE, to obtain appropriate subgrid scale fields. Thus this approach 
systematically implements model closures for macroscale discretisations. 

In Section El we prove that classic interpolation connects accurately the 
teeth across the gaps for the general linea r fourth order PDE. The technique 
of holistic discretisation ()R,obertjl2nn"il e.g.) resolves subgrid sca le struc- 
tures to reproduce with high fidelity the dy namics of specified pdes ( R.obertsI 



The techniques were adapted bv iRoberts &: KevrekidisI ()2005l ) to the 



gap-tooth scheme in the case where the microsimulator requires Neumann 
boundary conditions of specified slope/flux (see Sectional). Using the same 
techniques. Section analyses the use of classic interpolation of macroscale 
grid values in the microscale simulators and shows the following desirable 
properties: 

• the approach generates macroscopic discretisations which are consistent 
with the microscopic dynamics to high order in the macroscopic tooth 
separation H; 

• the macroscopic model and the microscopic solution field are essentially 
independent of the size of the teeth, measured by r; and 

• the macroscopic model and the microscopic solution field are essentially 
independent of the details of the tooth boundary conditions (tbcs) that 
couple the teeth together. 

Thus our proposed rule generates gap-tooth schemes that may be systemati- 
cally refined to high order accuracy, and gives rise to macroscale simulations 
that are largely independent of irrelevant microscale parameters. 



2 Dirichlet teeth (specified u) 

In this section we consider the case of microsimulators that require at each 
time step the field values on the edge of each spatial patch to be specified. 
We model this case by pdes with Dirichlet conditions coupling the dynamics 
in the teeth. The values for these Dirichlet conditions are obtained by inter- 
polation across the gaps between the teeth using finite difference operators 
and the exact relationships between the operators. When applied to simple 
diffusion, the resulting scheme has high order accuracy. 
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Discrete operators are essential in the analysis. Define the shift operator 
Eu(x) = u(x + H) and equivalently EU-j = Uj+i as appropriate for steps 
on the coa rse grid size H. Then we use the fo llowing identities for discrete 



operators (jNational Physical Laboratorvlll96ll . p. 65, e.g. 



mean yi = ^V^^ + E'^^^) , (1) 



difference 6 = E^/^ - E"^/^ , (2) 



shift E - 1 + ^5 + (3) 
derivative Hd^ = Zsinh"^ ^5 = 5 - ^5^ + 0(6^) , (4) 
yl^ = ^+lS\ (5) 

Formulae involving these operators become more accurate as the differences 5 
become small. Such small differences arise either as the macroscopic grid size 
H — > or equiva lently as the gradients of the physical field u become small. 
For example, [Roberts Sz Kevrekidid ^200Fh showed arbitrary order consis- 



tent macroscopic dynamics from a gap-tooth scheme as the grid size H ^ . 
The key to that analysis is the following transformation of the operator for 
evaluating spatial derivatives H9x at the patch boundaries E^^: 

F±^Ha^ = (1 + ^5 + ^5^)^^2sinh"^ ^6 . 

But this right-hand side, when expanded in a Taylor series in small differ- 
ences 6, is composed of terms which have an odd number of centred operators 
6 and yi. Consequently the right-hand side above would require field values 
halfway between the grid values. These are unknown. Instead, from 



multiply the right-hand side by the identity ii/a/I -|- 6^/4, and then expand 
in small differences 5: 



E±^Ha^ = (1 + ^6 + ^5^)±"2sinh"^ ^5 

|x5 + ^6^]^^2sinh-^ ^6 



^'■SO s' ' 24' J^*^ ' ^90 36' ^ 120 ' J" 

140 240 ' ^ 72' 720 ' J^" 
T^'^-SeO 1440 ^ 480 5040 ^ '^{'^ J ■ {^J) 

For microsimu lators with Dirichlet bounda rv conditions, we adapt the 
earlier analysis of R.oberts fc KevrekidisI ( 2nn5| ). But instead of determining 



the slopes at the tooth boundaries as above, the following interpolation of 
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Table 1: Growth rates A of perturbations from steady state u — 0: for 
diffusion (jH)) with m teeth, H = Irc/vn; with gap to tooth ratio r = 0.1 ; 
n — ^^ points in the microscale grid; and with the fourth order tbc 



m 


1 


2,3 


4,5 


6,7 


m + 1 : 2m 


4 


6- 10-12 


-0.946256 


-2.166285 


n/a 


-397.2 


8 


-3-10-12 


-0.996073 


-3.785024 


-7.121435 


-1588. 


16 


-1 - 10-1° 


-0.999750 


-3.984293 


-8.832102 


-6355. 


32 





-0.999986 


-3.998999 


-8.988613 


-25421 . 



the macroscopic field determines the field values u on the edges of the teeth: 

= (1+^6 + ^62^±r 

= 1 ir^tS + ^r^S^i lr(r2-1)^53 + lr2(r2-1)54 



± ^r(r2 - 1 ) - 4) + ^r2(r2 



5! ' 

± lr(r2 - 



1 



4)5^ 



1 



4)(r2-9)^5^ 
-4)(r2-9]6^ + 0(6^), 



(7) 



The pattern in the above interpolation formula is clear. Now we explore the 
numerical performance of a gap-tooth scheme using this formula to determine 
teeth boundary conditions. 

Consider gap-tooth simulations of the simple diffusion equation 



at 



a^u 
a^' 



and 27t-periodic in x. 



Imagine we only have access to the dynamics through a microscopic simulator 
of the diffusion (jHl), here coded by a fine discretisation on n grid points, spaced 
a distance r\ = H/(n — 1 ) apart, across a tooth of microscopic width h. = rH . 
The time integration is an explicit scheme with a microscopic time step, 
typically At = 10-^-10^. Fi gure 121 shows an example of the initially rapid 
microscopic evolution within one tooth; the microsimulator, coupled to its 
neighbors, rapidly evolves to a smooth state. Figure IHl similarly shows the 
initial evolution in two neighbouring teeth and how the smooth subgrid field 
arises through the coupling to the neighbouring teeth. Similar dynamics 
takes place during the initial instants of the Burgers' evolution shown in 
Figure [H 

Firstly we implement the following tbc. On the edge of the jth tooth, 
at X = Xj ± rH , the boundary condition of the fine discretisation is that the 
field 



Vj = [1 ± r[i8 + \r^8^ ± lr[r^ - 1 )^6^ + j^r^{r^ - 1 ]8^] Uj 



(9) 
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Figure 2: view of the initial microscopic evolution within a tooth with dy- 
namics described by the diffusion pde (jH)) and coupled to its neighbours. 




Figure 3: view of the initial microscopic evolution within a pair of neighbour- 
ing teeth with dynamics described by the diffusion pde (jHJ and also coupled 
to their neighbours. 
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Table 2: Growth rates A of perturbations from steady state u — 0: for 
diffusion (jH)) with m teeth, H = Irc/vn; with gap to tooth ratio r = 0.1 ; 
n — ^^ points in the microscale grid; and with the sixth order tbc. 



m 




1 


2,3 


4,5 


6,7 


m + 1 : 2m 


4 




5-10-'^ 


-0.981981 


-2.453767 


n/a 


-397.2 


8 


1 


•10-" 


-0.999653 


-3.927925 


-7.835158 


-1588. 


16 


8 


•10-" 


-1.000001 


-3.998611 


-8.966332 


-6355. 


32 


8 


•10-1° 


-1 .000002 


-4.000004 


-8.999518 


-25421 . 



The first few terms of d?)) provide this by interpolation from the surround- 
ing coarse grid values. For the jth tooth this tbc involves macroscopic grid 
values Uj_2, • • • , U-j+2 only, and thus we should be able to achieve 0(H'^) con- 
sistency with the microsimulator. We numerically linearize the map over one 
microscopic time step by systematically perturbing each and every micro- 
scopic value from zero (there are ran such microscopic values, one for each of 
n fine grid points in each of m teeth). We then transform the eigenvalues (x 
of this map to growth rates A = log((J.)/At. The mn growth rates fall into 
n groups of m modes. Each group corresponds to a microscopic internal 
mode of the dynamics; the mode is essentially the same in each tooth. Large 
negative growth rates correspond to rapidly decaying internal modes with 
significant microscopic structure within each tooth. The group of m modes 
with small growth rates correspond to the relatively slowly evolving macro- 
scopic modes of interest that arise through the coupling of the microscopic 
dynamics across the teeth. Table ^ shows the leading seven growth rates, 
and the magnitude of the leading internal growth rate, for various numbers 
of teeth, m 4, 8, 1 6, 32 . The exact growth rates of the diffusion pde (jH)) are 
A = — for integer k. The table shows that as the number of teeth doubles, 
the accuracy of the growth rates of the macroscopic modes improves by a 
factor of about 16. This is consistent with an C(H^) method as predicted 
for diffusion with tbc (0). 

Table |21 shows the even higher order accuracy from implementing sixth 
order TBCs from (0) — growth rates slightly larger than the ideal seem to be 
due to the relatively small number of microscopic grid points within the teeth. 
These sixth order tbcs are used in the simulations of the nonlinear Burgers' 
equation shown in Figure^ This simulation suggests that gap-tooth schemes 
employing such tbcs even for nonlinear systems are effective. 
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Table 3: Growth rates A of perturbations from steady state u — 0: for 
diffusion (jH)) with m teeth, H = In/m.; with gap to tooth ratio r — OA ; 
n = 1 1 points in the microscale grid; and with the mixed tbc fll()|l with 



a = 0.95 and b ^ 0.05 . 



m 




1 


2,3 


4,5 


6,7 


m + 1 : Im 


4 


9 


•10-1^ 


-0.939448 


-2.151772 


n/a 


-240.7 


8 




2- 10-" 


-0.990854 


-3.766007 


-7.089004 


-756.9 


16 


7 


•10-" 


-0.996405 


-3.971027 


-8.803476 


-2417. 


32 


6 


•10-1° 


-0.998047 


-3.991243 


-8.971213 


-8153. 



3 Mixed boundary conditions for the teeth 

Let us explore mixed boundary conditions at the edges of the teeth: suppose 
the microsimulator requires avj ± b9xVj specified on the edge of the teeth 
X = Xj ± rH for some constants a and b. The case a = 1 and b = 
constitutes Dirichlet tbcs discussed in the previous s ection. The case a = 
and b = 1 constitutes Neumann tbcs as discussed by Roberts &: KevrekidisI 



(200i): there we used the interpolation formula (0) to specify slopes/fluxes 



on the edge of each tooth; we obtained spectra of accuracy similar to those 
in Tables ID and 

For mixed tbcs we propose to simply combine (jH)) and (|7j) to give, for 
example, the fourth order in macroscopic grid size H = In/m boundary 
condition 

avj ± bQxVj 

= a [1 ± r^5 + {r^8^ ± lr{r^ - 1 )^5^ + ^r\r^ - 1 ]b^] Uj 

± ^ [^t5 ± rb' - (I - lr^)yi8' T - |r^)6^] Uj 
on X = Xj ± rH . (10) 

We use a = 0.95 and b = 0.05 in the mixed tbc: this gives a mixed boundary 
condition where the effects of the function value Vj and its gradient d^Vj are 
roughly comparable in the TBC (if the parameter b is significantly larger, then 
the gradient term dominates the tbc). The numerical eigenvalues given in 
Table Elfor the diffusion equation (jHJ with these tbcs again show convergence 
to the correct eigenvalues as the number of teeth increases, that is, as the 
macroscopic grid size H ^ . However, the convergence is not as rapid as 
for Dirichlet tbcs. The poorer convergence as H — > seems to be due to 
the microscopic grid resolution: successive doubling of the number of interior 
points, see Table EJ demonstrates that there are significant errors of C'('n^) 
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Table 4: Growth rates A of perturbations from steady state \i — 0: for 
diffusion (jH)) with m = 8 teeth; with gap to tooth ratio r = 0.1 ; n points 
in the microscale grid to show variation with microscale resolution; and with 
the mixed tbc dTUl) with a = 0.95 and b = 0.05 . 



n 


1 


2,3 


4,5 


6,7 


11 


-2- 10-" 


-0.990854 


-3.766007 


-i.omoA 


21 


1 • 10-1° 


-0.99AS9G 


-3.781289 


-7.117588 


41 


-3-10-1° 


-0.995792 


-3.784677 


-7.123924 



in the microscale grid size r|. Thus the total error in this implementation of 
the mixed tbcs seems to be 0(H'^,ri^). 

Here the microscale simulation is that of a fine discretisation of a pde. 
Thus the derivatives in the mixed tbc (fTn|) are subject to the significant 
errors of numerical differentiation when computed on the microscale. As Ta- 
ble m shows, the approximation of derivatives does incur errors; we would 
be better off without such errors. Higher order formulae for microscale in- 
terpolation would reduce the microscale errors in the boundary derivatives, 
perhaps from C(ri^) to 0{r\^^^ but would ruin the small bandwidth of the 
microscale simulator. In any case, recall that we adopt the policy that we 
cannot change the microscale simulator as it is a legacy code handed to us 
from past development. We cannot (do not want to) change the nature nor 
accuracy of its boundary conditions. Consequently we proceed to address the 
problem of supplying boundary conditions at the edge of the teeth, precisely 
as required during execution of the legacy microscale simulator. 

4 Teeth with two point boundary conditions 

A microscale simulator may have implemented boundary conditions that do 
not fit into the classic partial differential equation form of Dirichlet, Neumann 
nor mixed. Here the microscale simulator implements a discretisation of 
3 point stencil width. Consequently the simulator has been written so that 
the supplied boundary conditions only depend upon each of the two extreme 
pairs of points in each tooth. We thus investigate teeth boundary conditions 
that specify a combination of these two point values of the field at the edge 
of each tooth. This specific case is just one example of the wide range of 
possible nonlocal tbcs that specific microsimulators may require. 

Suppose the microsimulator, here a fine spatial discretisation of the dif- 
fusion PDE (jHl), implements a tooth boundary condition of the (linear) form 
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Table 5: Growth rates A of perturbations from steady state u = : for 
diffusion (jH)) with m = 8 teeth; with gap to tooth ratio r — OA ; n points in 
the microscale grid to show variation with microscale resolution r\ oc 1/n; 
and with the fourth order general tbc ()12|) with (3 = 1. 



n 


1 


2,3 


4,5 


6,7 


11 


1 • 10-'" 


-0.999741 


-3.984137 


-8.831209 


21 


-3-10-1° 


-0.999742 


-3.984159 


-8.831368 


41 


1 • 10-^ 


-0.999742 


-3.984169 


-8.831453 



Vj 1 + (3vj 2 and (3vj n-i +Vj,n are specified, (11) 

where Vj_i denotes the microscale field value at the ith microscale grid point in 
the jth tooth. For example, the case (3 = 1 approximates Dirichlet boundary 
conditions at the microscale grid mid-points and Xj^n-i/z near the edges 
of each tooth, which is exactly the case a = and b = 1 implemented in the 
previous Sectional Different values of (3 would approximate different mixed 
boundary conditions of the previous section. 

The procedure is straightforward: we interpolate the macroscale grid val- 
ues to find the specific values required by the boundary conditions (fTT|) . 
Recall that gives a fourth order interpolation from the macroscale grid 
to points at the tooth boundaries x = Xj ± rH ; this gives appropriate values 
for Vj 1 and Vj^ri- Get appropriate values for Vj 2 and Vj^n-i through sim- 
ply replacing in the formula the ratio r = h,/(2H) by the ratio required to 
reach the penultimate microgrid point, namely r' = (h,/2 — ri)/H, where 
r| = h,/(rL — 1] is the microgrid size. Thus the fourth order version of the 
boundary condition (jlip is that at x = Xj ± rH 

(1 + (3E^^/")vj 

= [1 ± r^5 + lr^8^ ± ^r(r^ - 1 + ^r^(r^ - 1 jS^] Uj (12) 



1 ± r'yib + Ir'V ± ^r'(r'^ - 1 ]^6^ + ^r'^(r'^ - 1 ]6 



Implementing the TBC p2|l for the diffusion equation (jH)) gives a numerical 
approximation scheme with eigenvalues shown in Table El for varying mi- 
crogrid resolution. See that there is only an extremely weak dependence 
upon the microgrid size r\. Thus implementing directly the boundary condi- 
tions that the microscale simulator actually expects during execution results 
in much better accuracy than trying to approximate the microscale TBCs us- 
ing computed spatial derivatives. 
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Table 6: Growth rates A of perturbations from steady state u — 0: for 
diffusion (jH)) with m teeth, H = Irc/vn; with gap to tooth ratio r = 0.1 ; 
n — ^^ points in the microscale grid; and with the fourth order general 
TBC (fnil with |3 ^ 1 ■ 



m 




1 


2,3 


4,5 


6,7 


m + 1 : 2m 


4 




8-10-'^ 


-0.946069 


-2.165068 


n/a 


-489.5 


8 


4 


•10-" 


-0.996034 


-3.784277 


-7.118312 


-1958. 


16 


1 


•10-1° 


-0.999741 


-3.984137 


-8.831209 


-7832. 


32 


8 


•10-1° 


-0.999983 


-3.998964 


-8.988427 


-31329. 



Lastly, Table El shows the eigenvalues of the gap-tooth scheme for varying 
number m of teeth in the domain. See that the eigenvalues converge to their 
correct values like C(H'*) as expected by the construction. 

Higher order tbcs, in the macroscopic grid size H, would similarly be 
based upon the expansion We then expect even more rapid convergence 
as the macroscale grid size H — > . 

5 The model is independent of the tooth bound- 
ary conditions 

Here we use analytic methods of holistic discretisation 

to explore the gap-tooth scheme on a general class of pdes with general mixed 
boundary conditions. The analysis establishes three important properties: 

• the approach generates macroscopic models which are consistent with 
the microscopic dynamics to high orders in grid spacing H; 

• the macroscopic model and the microscopic solution field are essentially 
independent of the size of the teeth, as parametrised by r; and 

• the macroscopic model and the microscopic solution field are essentially 
independent of the details of the TBCs. 

5.1 Theory underpins analysis of a PDE with tooth 
boundary conditions 

We explore solutions of the class of linear hyper-advection-diffusion pdes 
9u 9^u / 9u , 9^u 9%\ 



RobertsI jml e.g. 
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where a, b and c are arbitrary parameters, and where e is introduced as a 
convenient mechanism to control truncation in the multivariate power series 
solutions in the parameters a, b and c. This pde is solved with mixed tooth 
boundary conditions inspired by (fTn|) . namely that on x = Xj ± rH, and in 
terms of an artificial parameter y that we explain shortly, 

= ±cx {1 + yr [±^5 + {rS^] + y^r[r^ - 1 ) [±^5^ + ^rS^] 

+ 1 {y [^15 ±r6^] [-(I - lr')^8'Tr[^, - -^jS"] 

+ y'H-y + ^/)^6^ ± - + ^r^)5T } 
+0(y4). (14) 

Explore the structure of this complicated looking tbc: iavj + S^Vj represents 
a general linear combination of the microscopic field at the edge of each tooth 
that needs to be specified for the microscopic simulator; those terms in the 
right-hand side multiplied by ±(x form the estimate of the field Vj interpolated 
from the surrounding macroscopic grid values; those terms in the right-hand 
side multiplied by 1 /H form the estimate of the field's gradient 9xVj inter- 
polated from the surrounding macroscopic grid values. However, these two 
interpolations only hold when the artificial parameter y = 1 ; one is the phys- 
ically interesting value of y. Why then do we introduce the parameter y? 



The reason is that, as in "discretisation" ()Roberts.2001ai) . based around the 
special values of the parameters y — (x — e = , the general pde fjl3|) with 
TBC ()14|1 possesses a (slow) centre manifold parametrised by the macroscopic 
grid values. On this centre manifold the evolution of these macroscopic grid 
values forms a macroscale model of the pde. This model has rigorous theo- 
retical support based upon y — , and it becomes physically relevant when 
evaluated at y = 1 . 

We briefiy explain how centre manifold theory underpins the macroscale 
model. Initially set y = a = e = 0; then the pde+tbc become the diffusion 
equation with insulating boundaries at the edges of the teeth, x = Xj ± rH . 
Thus, exponentially quickly, all structure within each tooth diffuses away to 
become constant, but a different constant for each tooth depending upon 
the initial conditions. See a similar evolution in Figures |21 and El but there 
the teeth are coupled, so that the rapid evolution is to a smooth variation 
in each tooth, whereas here the insulated evolution, y = a = e = 0,isto 
a constant in each tooth. But we are only interested in fully coupled teeth 
for which y = 1 , and in non-zero a and e. Thus from the simple base of 
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piecewise constant fields, we construct a description of the field u and its slow 
evolution as a power series in the "perturbations" measured by y, <x and e. 
The departure of the field u from a constant within each tooth gives the 
microscopic (subgrid, subtooth) field, as shown for example in the smooth 
fields of Figures El and El that are quickly established. The slow evolution of 
the coarse grid values Uj gives the macroscopic model. 

The various powers of y in the tbc ()14j) are chosen so that trunca- 
tion of the expressions to errors 0(y^) will generate a discrete macroscopic 
model expressing Uj in terms of only Ui_ p.^i , . . . , Ui+p-i (a spatial sten cil of 
width 2p — 1). Centre manifold theory ( Can 198ll Kuznetsov 19951 e.g.) 
asserts that 

• such a model exists, 



• that through its exponential attractiveness, the model is relevant in 
some finite neighbourhood ofy = (X=e = 0, 

• and that we may systematically construct the power series approxima- 
tion to the model. 

Because truncation to errors O(y^) results in a model with stencil width 2p — 
1 , such a truncation corresponds to the gap-tooth scheme utilising tbcs in- 
volving interpolation from only the 2p— 1 neighbouring grid values Uj_p+i , . . . , Uj+p_i . 
Cor nputer algebra^ performs all the tedious details of constructing the 



model ()Roberts.l997.) . We seek a model where the subtooth/subgrid field 

u(x,t) = Vj(x,ll; a, y, e] . (15) 

That is, the subtooth field has some spatial structure, such as that in Figures 
ElandEl which: depends upon the neighbouring grid values Uj_p+i , . . . , Uj+p_i ; 
may depend upon the specific tbc through its parameter a; depends upon 
the specific pde though its parameter e; and depends upon the coupling 
parameter y. Centre manifold theory assures us the evolution of the system 
is governed by the evolution of the grid values: 

Uj = gj(U;cx,y, e); (16) 

this formula is the macroscopic (closed) discretisation. We solve by iteration 
the PDE (Uni) with TBC (dH) to find the centre manifold (fT3|) and its associated 
coarse discretisation ()16p. The results are expressions for the microscopic 
fields Vj and the macroscopic evolution, Uj = g-j , that are accurate to some 
specified order in the small parameters a, y and e. 

^http : //www . sci .usq. edu. au/staf f /aro ber'ts/CA/burgennixed. red| is the source 
script which was available at the time of writing. 
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5.2 Modelling 0(e) changes to the PDE 

The subgrid fields in each tooth will in general depend upon the coefficients 
a, b and c that determine the PDE. For example, there are nontrivial depen- 
dencies upon the advection speed that ensure the macroscale model naturally 
trans forms to an upwind discretisation for large advection speeds c (Roberta 



2002[ ) — such influences show up in the C'(e^) terms that we explore in the 



next subsection. Here we ffist explore the models linear in a, b and c, that 
is, linear in the general changes to the PDE ()13|1 with the tbc (|T^ . 

For example, to errors 0((X^,y'*, e^), computer algebra generates the 
macroscopic evolution 

When evaluated at the physically relevant parameter y = 1 these are the 
classical finite difference operators for the PDE (fT^ . truncated to C(5^). 
Consequently the terms in the macroscopic model ()17p are consistent with 
the PDE p3|) to various orders in the macroscopic grid size H. The order of 
consistency depends upon the order of truncation in the artificial coupling 
parameter y and the order of the derivatives in each term. Observe that the 
macroscopic evolution operator is independent ofr, the size of the teeth, and 
independent of oc which parametrises the precise nature of the tbc ^14\ )- 

Now we explore the microscopic field within the teeth. To low order in the 
coupling parameter y and in terms of the microscopic tooth space variable 
£,= (x-Xj)/H, we find 

= U^+y{[£,^5 + l£.V] 

+ ecH [(^£,3 - {r^Q + iHocr^E^ - ^H^cxVE,] 5^} Uj 

+ C(a^y^e2). (18) 

The ffist line gives the classic quadratic interpolation through the grid val- 
ues Uj and Uj-i-i. The second line shows microscopic field structure in the 
advection speed c. But it exhibits undesirable dependence upon the tooth 
width r and nature a of the tbc. However, inspect the next order terms in 
coupling parameter y: 
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+ C(a^y^e2); (19) 

The dots denote the terms given in the right-hand side of (fTHj) . The first 
hne in the above higher order terms contains reassuringly the classic quartic 
interpolation formulae. The second line, when we set y = 1 , cancels all 
the undesirable a and r dependence in the lower order (fH?]) . The third 
and later lines above describe higher order microscopic structure; some of 
this undesirably depends upon the tbc through a and the size of the teeth 
through r, but as far as we have explored, any a and r dependence introduced 
at any order in y is canceled by terms at higher orders in y. Thus any finite 
truncation of the power series expansion for the model may have undesirable 
dependence upon a, and r, but as the order in the artificial parameter y is 
increased, this dependence is removed. In this sense, the microscopic field is 
"essentially" independent of the details of the tbc and independent of the 
tooth width r. 

5.3 Modelling 0{e^) effects in the PDE 

In our analysis we find that even order operators in the microscale PDE, 
such as the diffusion u^x and the hyper-diffusion au^xxx, are represented 
simply in the macroscale discretisation. However, odd order operators in 
the microscale PDE, such as advection cUx and the dispersion buxxx, create 
nontrivial effects; these first show up in terms quadratic in their amplitude 
and hence they first appear in terms of O(e^). We now show that the 
effects of C'(e^) terms typically act to stabilise the C(e) discrete model. The 
implication for the gap-tooth method is that the resolution of the subgrid 
structures by the microscale simulator will also typically maintain stability of 
the discrete macroscale model. That is, the microscale simulator will provide 
successful closure for the macroscale discretisation when coupled with the 
proposed tbcs. 

To illustrate this, we can construct the approximate model of the PDE (fT^ 
with TBC (HD) to errors C(cx^,y'^, e^); that is, we include quadratic effects in 
the coefficients a, b and c. The details of the model are too long to record 
here. However, we find that the equivalent PDE to the macroscale discrete 
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model ()17p with its C'(e^) modifications is 



9u 9^u 



9t 



9x2 



9u 7^ 9^u 
'YC—+y b— - 
9x ox^ 



94u 



9x4 



-(y - y^)bc— - - (y^ - y^jb^ 



+ H'|i^(y-y 



9x2 
9% 



9x4 



+ e 



4(y-y')c 



9x4 



9x3 



(y'-y']b 



9^u 
9x5 



^(y'-y')a 



9^u 
9^ 



4(y-y')cV 



9x2 



ll(y-5y^ + 4y3)bc 



,94u 



+ |(y' - T')bcr^ - -^(y^ - 7y3 + 4y\^)b^ 



9x4 

9^u" 
9x6 



(20) 



Observe that to this level of accuracy there is no dependence upon the tbc 
parameter a, thus our comments apply for all the tbcs. Now consider the 
components of ()2()|1 in turn. The first line of (pn|) is the original general linear 
PDE (fT!^ when evaluated at the physically meaningful y = 1 . The second 
line shows some error terms, quadratic in e, that disappear for y = 1 . The 
bcUxx error disappears when O(y^) terms are retained, which is as soon as 
the discretisation stencil is wide enough to model the third order dispersion 
term buxxx- The b^Uxxxx dispersion-induced term shows that the method 
initially incorporates its effects as enhanced dissipation, as the coefficient 
of y^ is negative; then, when higher order accuracy is requested by retaining 
0(y3) terms, the method proceeds to remove the incurred error via the 
y3 term. 

Consider applying these tbc to a microscale simulator; for slow enough 
spatial variations, the microsimulator is equivalent to some 'infinite order' 
PDE. For example, the microscale discretisation — (1/ri^]6^Ui is, by 
equivalent to the PDE Ut = (4/ri^] sinh^(r|9x/2)u. We expect the behaviour 
of errors in the gap-tooth scheme seen here, induced by the Uxxx and Uxxxx 
terms, be representative of the behaviour of errors in the 'high order' equiv- 
alent terms of any given microscopic simulator. 

The remaining terms in ()2U|) are C(H^) and hence vanish as the macro- 
scopic grid size H — > . Nonetheless look at the terms in the third to sixth 
lines as they apply to simulations with finite H. The third and fourth lines 
show that the initial discretisation errors of the linear terms are eliminated, 
for the physical y = 1 , via the next higher order in coupling parameter y. 
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The particular focus of this subsection is the terms on the next two hues. 
The c^Uxx terms show that at low order truncations in y the method treats 
advection in a manner that increases dissipation, as the coefficient is posi- 
tive, and thus helps to maintain the sta bility of the macroscale discretisation 
to high advection speeds c (explored in lRobertsll2n(mll2nn2h . The interac- 
tion between advection cUx and dispersion bu^xx can maintain stability or 
be destabilising depending upon the sign of be: whether we truncate at y 
when — bcUxx dominates (second line) or we truncate at y^ when bcUxxxx 
dominates (fifth and sixth lines), the combination is stabilising whenever 
be < , that is, when the phase velocity of wave-like effects does not change 
direction as a function of wavenumber. The last term on the sixth line will 
be dominated by the dissipative b^Uxxxx term on the second line, and we 
presume will vanish at higher orders in the coupling parameter y. Thus, 
from the equivalent pde ()2()j) of the macroscale model, we deduce that the 
subgrid scale interactions between processes in the pde, and hence for mi- 
croscale simulators in general, are accounted for in this approach to generate 
a macroscale model that is typically stable. 

Indeed this equivalent pde ()20j) confirms support for the gap-tooth scheme 
with TBC by centre manifold theory. Theory asserts that the original system, 
here the pde (jSDI), and the centre manifold model, here the macroscale dis- 
cretisation (fTH|l . have the same stability. Thus when the microscale system is 
stable, so will the macroscale discretization. The caveat is that we can only 
construct the centre manifold approximately; we control the errors to some 
order in the parameters a. y and e, but there will be some error, albeit of 
high order in the parameters. 

6 Conclusion 

We use macroscale interpolation based upon the expansions (jHl) and ((Zj) to 
determine tbcs for the boundary conditions at the edge of the teeth in the 
gap-tooth scheme. The interpolation was used to implement directly what- 
ever boundary conditions are actually needed by the microscale legacy code 
during execution. Figure |3] shows a simulation of the nonlinear Burgers' 
equation with 2 point boundary conditions at the boundaries of each tooth 
as an illustrative example. We found that the macroscopic models resulting 
from the microsimulator and the constructed TBC were consistent, to high 
order, with the microscopic dynamics; that the macroscopic models and the 
microscopic (subgrid, subtooth) fields were essentially independent of the 
tooth size and the detailed nature of the TBC. We expect the same type of 
TBC to be effective for microsimulations in more than one spatial dimension. 
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Figure 4: simulation of Burgers' equation using general 2-point boundary 
condition on the teeth of the fourth order ()12|1 with |3 = 1 demonstrates the 
method is stable even for nonlinear pdes. 

Interesting future research would seek tbc that do not require communica- 
tion across the gaps between the teeth at each and every microscale time 
step, and the interplay of TBCwith implicit integration schemes. 

Further exciting research would explore issues of existence and perfor- 
mance of TBCs for stochastic microsimulators. 
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